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ABSTRACT 


A  new  approximation  technique  to  a  certain  class  of  nonlinear  filtering 
(signal  processing)  problems  is  considered  here.  The  method  is  based  on  an 
approximation  of  a  nonlinear,  partially  observable  system  by  a  bilinear 
model  with  fully  observable  states.-  The  filter  development  proceeds  from 
the  assumption  that  the  unobservable  states  are  conditionally  Gaussian  with 
respect  to  the  observation  initially.  The  method  is  shown  to  be  promising 
for  real-time  communication  and  sonar  applications  as  demonstrated  by  com¬ 
puter  simulations.  Moreover,  some  of  the  traditional  techniques  evolve  as 
special  cases  of  this  methodology. 
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I.  INTRODUCTION 


In  general,  the  optimal  nonlinear  mean  square  filtering  (signal  process¬ 
ing)  problem  does  not  have  a  finite-dimensional  recursive  synthesis.  This 
is  due  to  the  fact  that  the  conditional  probability  density  of  the  unob¬ 
servable  states  with  respect  to  the  observations  cannot,  in  general,  be 
characterized  by  a  finite  parameter  set,  such  as  for  a  Gaussian  distribution 
utilizing  linear  methods.  Thus,  approximation  and  ad-hoc  techniques  must  be 
employed  to  construct  practical  filters  for  nonlinear  systems. 

In  this  paper  a  new  nonlinear  finite-dimensional  filtering  methodology 
is  presented.  The  method  which  could  be  called  a  two-step  nonlinear  filter 
(TNF)  is  of  a  more  general  nature  than  the  linear  model  approach  (or  extended 
linear  filtering)  ,  and  does  not  require  the  model  "smoothness"  which  is 
crucial  to  most  of  the  existing  techniques.  Hence,  the  proposed  technique 
may  expand  the  range  of  practical  problems  that  can  be  handled  by  nonlinear 
filtering. 

Herein,  the  general  nonlinear  model  is  approximated  by  a  control  model 
of  "bilinear  form."  The  "best"  model  approximation  of  this  form  is  then 
computed  by  the  appropriate  stochastic  control  [1].  The  final  step,  with 
the  feedback  control  as  a  function  of  the  estimated  state  x^,  requires  com¬ 
putation  of  the  optimal  (m.  s.e.)  state  estimator  by  the  conditionally 
Gaussian  filter  which  is  formally  developed  by  Liptser  and  Shiryayev  [2] 
and  extended  by  Kolodziej  [3]. 

Extension  of  the  results  to  large-scale  nonlinear  systems  is  accom¬ 
plished  by  incorporating  a  novel  decomposition  scheme  in  the  filter  design 
to  alleviate  the  complexity  of  the  control  problem  and  the  "curse  of 
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dimensionality"  of  the  filtering  algorithm.  Formal  representation  of  the 
filter  is  given  in  Figure  1. 

Application  of  the  developed  filter  to  a  scalar  nonlinear  system  which 
lacks  model  smoothness  is  discussed  and  application  of  the  derived  multi¬ 
dimensional  filtering  algorithm  to  a  low-order  nonlinear  tracking  problem 
according  to  a  global  criterion  is  presented.  In  addition,  a  comparison 
with  traditional  methods,  such  as  the  popular  Extended  Kalman  Filter  (EKF) 
is  given  by  digital  computer  simulation  to  demonstrate  the  effectiveness 
of  the  obtained  results. 
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The  dashed  tines  implies  functional  representation. 
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II.  MODEL  APPROXIMATION  AND  FILTER  STRUCTURE 

Consider,  as  given,  some  complete  probability  space  (ft,F,P).  Let  (w^) , 
i  =  1,2  defined  on  (ft,F,P)  be  mutually  independent  Wiener  processes.  Also, 
given  the  following  nonlinear  models,  which  may  represent  a  broad  class  of 
nonlinear  systems,  with  dynamics  described  by  a  family  of  stochastic  differ¬ 
ential  equations  (in  the  Ito  sense)  of  the  form 

dxt  -  F(x,y, t)dt  +  G(x,y,t)dw^  , 

(1) 

2 

dyt  =  H(x,y, t)dt  +  R(x,y,t)dwt  , 

where  (F(*)>  H(«))  are  nonlinear,  real  vector  functions,  and  (G ( • ) ,  R(*)) 
are  matrices  of  compatible  orders. 

The  optimal  m.s.e.  filter  for  the  system  in  (1)  is  known  to  be  the  con¬ 
ditional  expectation  of  the  system  given  the  observation  (y  ;  0  s  t), 

s 

t  e  [0,T]  provided  that  a  finite  second  moment  solution  exists  with 

xt  =  E(xt/yg;  s  e  [0,t])  .  (2) 

In  principle,  a  sequential  version  of  (2)  can  be  found,  but  in  general, 
the  recursive  formulae  consist  of  an  infinite-dimensional  system  of  moment 
equations  which  are  needed  to  characterize  completely  the  conditional  proba¬ 
bility  density,  p(x  ,t/y  ;  s  £  [0,t]).  Thus,  one  is  forced  to  seek  an 

t  s 

approximation  technique  for  practical  implementation. 
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Herein,  a  new  technique,  which  generates  a  finite-dimensional,  nonlinear 
filter  and  a  "close"  (m. s.e.)  model  approximation  to  the  original  model  in 
(1)  is  presented.  This  technique  may  be  called  "an  approximation  in  the 
model  parameter  space,"  and  its  parameters  are  functions  of  the  feedback 
control  law,  u^,  which  is  itself  a  function  of  the  observation  process  y  . 

In  general,  the  model  in  (1)  can  be  approximated  as  follows: 


dx  =  F (x  , u  , t)dt  +  G(u  ,t)dw1  , 

L  L  U  L  t 

dy  =  H(x  ,u  , t)dt  +  R(u  ,t)dw2  , 

*-  Cl  t  t 


(3) 


where  (F,G,H,R)  are  of  appropriate  dimensions  and  are  functionals  of  u  . 
The  control  u^  is  a  measurable  stochastic  process  with  respect  to  Y  ~ 
a-algebraCy^;  s  e  [0,t],  and  it  is  chosen  to  minimize  the  following  mean- 
square  (global)  criterion: 


Q(u)  =  E (  /  ||k  -  k  1  2  dt)  , 

0  C  t 


(4) 


where  kt  denotes  any  of  the  functions  F,  H,  R,  or  G  in  (1),  and  kt  is  its 
corresponding  approximation  in  (3) .  Here  ||  •  ||  is  the  Euclidean  norm,  and 
the  arguments  (x^y^jt),  (x^jU^jt)  are  omitted  for  brevity.  It  is  note¬ 
worthy  that  the  choice  of  (4)  ensures  a  form  of  global  filtering  criterion 
even  though  the  search  for  the  approximation  k£  falls  into  a  class  of  sto¬ 
chastic  control  problems  and  depends  strongly  on  the  type  of  nonlinearities 
in  the  system. 


An  important  special  class  of  (3)  is  the  following  bilinear  form  ; 
dxt  -  [A(ut,t)xt  +  B(ut,t)]dt  +  G(ut,t)dwJ  , 


dyt  ~  [C(ut’t)xt  +  D(ut,t)]dt  +  R(ut,t)dwJ 


(5) 


Where  (A,B,C,D,G,R)  are  linear  fnnctionals  of  ^  and  are  of  appropriate  di- 
mensions.  The  term  bilinear  refers  to  the  fact  that  the  system  is  linear 
in  control  and  state,  but  not  Jointly  linear  (4).  It  has  been  shown  that 
such  bilinear  systems  are  quite  common  in  physical  applications  and  repre- 
sent  a  significant  class  in  their  own  right. 

Now  assume  v  which  minimises  (4),  is  available;  then  (5)  is  a  "close 
“  the  sense  of  (u) ,  approximation  model  to  system  (1)  since  the  minimisa¬ 
tion  criterion  is  a  measure  of  the  quality  of  the  approximation.  If  it  is 

assumed  that  the  initial  unobservable  state  x  r*r*n^*+-* 

9  Xq,  xs  conditionally  Gaussian 

with  respect  to  the  initial  observable  state,  yQ,  then  under  certain  broad 

assumptions  (see  [2,  or  (3])  the  optimal  ......  filter  for  the  system  in 

(■5)  takes  the  form  of 


dxt  =  (Axfc  +  B)dt  +  Advt  , 

drt  =  (AFt  +  FtA*  +  G  G*  -  AA*)dt  , 
A  =  (I*  C*)!"1  , 

dvt  =  *  1tdyt  -  (Cxt  +  D)dt]  , 

-2  '  - 
R  =  R  R*  , 

where  x  -  E(xt/Yt),  ^  .  Cov  (x^)  . 


(6) 
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Here  E(/)  denotes  the  conditional  expectation  operator,  T  conditional  covari¬ 
ance,  *  transposition  operator,  and  is  the  a-algebra  generated  by  the 
observations  on  [0,t],  t  e  [0,T].  Again  the  arguments  (u^jt)  are  omitted 
for  brevity. 

A 

Now,  if  the  optimal  control  or  its  approximation  can  be  obtained, 
then  system  (6)  provides  a  finite-dimensional  approximation  for  the  nonlinear 
filtering  of  (1).  Also,  the  partially  observable  system  in  (5)  is  trans¬ 
formed  into  a  completely  observable  system  (6).  Consequently,  to  solve  for 
the  control  law  analytically,  the  minimization  criterion  (A)  must  be  trans¬ 
formed  accordingly.  Thus,  if  (4)  is  rewritten  in  the  form 

T 

Q (u)  =  E  /  L (x  ,u  , t)dt  , 

0  c  C 

then  for  u  -  u 

X 

Q (u)  =  E(  /  L(x  ,T  ,u  ,t)dt)  ,  (7) 

0  etc 

^  ^  00 

where  L(x  ,T  ,u  ,t)  =  /  L(£,u  ,t)f(£,x,r  )d£,  and  f(£,x,T  )  denotes  a  dif- 

t  t  t  t  t  t 

ferential  Gaussian  measure  with  mean  x^,  variance 

If  u  is  replaced  in  (6)  by  u,  then  an  equivalent,  completely  observable, 
stochastic  control  problem  [(6),  (7)]  emerges  where  the  new  state  of  the  sys¬ 
tem  (x^jT^)  are  generated  by  (6).  Hence,  the  filtering  problem  is  actually  re¬ 
placed  by  a  stochastic  control  problem  which  results  in  a  stochastic  non¬ 
linear  Bellman  equation  [1].  Apparently,  this  could  lead  to  a  more  difficult 
problem  to  solve  than  the  original  filtering  problem  for  (1).  However, 
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satisfactory  approximation  to  the  control  law  can  be  found  without  solving 
the  Bellman  equation  exactly. 

It  might  be  noted  here  that  the  solution  to  this  problem  according  to 
ut  resembles  the  approximation  criterion  in  the  EKF  approach,  and  makes  the 
approximation  of  closed  analytical  form. 

The  following  scalar  example  illustrates  the  suggested  filtering  proce¬ 
dures. 


Example  1 

Consider  a  process  with  an  absolute  value  detector  so  that 


dx,  = 


dy+.  = 


f |xt|dt  +  a1  dwfc  , 


hx  db  +  a„  dw  , 
t  2  t  ’ 


(8) 


1  2 


where  wt>  are  independent  Wiener  processes,  f,h,a^,a2  are  constant,  and 
t  e  [0,T].  The  bilinear  feedback  approximation  yields 


dx  -  u  x  dt  +  a,  dw  , 
t  t  t  It* 


dyt  =  xtdt  +  a2  dwt  , 


(9) 


where  u  is  selected  to  minimize 


E(  f  (f | x  |  -  ux  )  dt) 
0  C 


(10) 
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The  optimal  (m.s.e.)  filter  for  (8)  with  xQ  conditionally  Gaussian  w.r.t., 
Yq  has  the  form  of 


-1 

dxt  =  UX1  dt  +  rt  h(a2^  dvt  * 

drt  =  t2Vt +  4  - h2  rt(c2)-2^dt  * 


(id 


where  dvfc  -  (a^  ) (dy^  -  hxt>dt  is  the  innovation  process.  Now  from  (7)  the 
equivalent  minimization  criterion  to  (10)  is 


Q(t,x,t)  =  E(  /  (u | C I  -  uC)2d$(£)  , 

0 

*  2 

where  d$(5)  =  ~  exp(-0.5  — )d£,  for  nonsingular  r  . 

.  ZTt  1 

found  that  the  optimal  u  is  given  by 


(12) 

From  [5]  it  is 


2xT  "2 

u  ~  f  [ - exp  (-  “■)- 

/2Trrt  2F 


2(x2  +  r)erf(-  — )  -  (r  +  0.5x  (x2  +  r)  1 

Sir  3r  3X 


(13) 


where  v  is  the  solution  of  the  stochastic  Bellman  equation  which  means  that 
we  have  to  solve  a  nonlinear  partial  differential  equation  first  in  order 
to  evaluate  the  optimal  control  in  (13).  However,  it  is  seen  that  an  approxi¬ 
mating  optimal  control 


u  =  f 


[2x1 exp (-x/2r) ] 

^2rrT  (x2+r) 


+  2  erf(-x//2 F) 


(14) 


already  gives  better  performance  than  the  EKF  [5]. 
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III.  DECOMPOSITION  SCHEME  FOR  MULTIDIMENSIONAL  SYSTEM 


Extension  of  the  previous  results  to  large-scale  nonlinear  systems  is 
accomplished  by  incorporating  a  novel  decomposition  scheme  in  the  filter 
design.  This  should  alleviate  the  "curse  of  dimensionality"  which  is  en¬ 
countered  if  the  optimal  approximation  is  applied  directly.  Recall  that  this 
solution  involves  a  function  of  (n  +  variables,  where  n  is  the  dimen¬ 

sion  of  the  system. 

The  strategy  adopted  in  this  decomposition  scheme  is  based  on  the  decom¬ 
position  of  the  original  system  into  two  interconnected  subsystem  approxima¬ 
tions.  The  first  subsystem  is  only  linear  (i.e.,  linear  process  and  observa¬ 
tion  equations).  The  second  subsystem  includes  all  the  nonlinearities  in 
the  system  which  are  then  approximated  by  the  proposed  TNF.  The  main  ad¬ 
vantage  of  using  this  scheme  is  that  the  control  parameters,  which  are  needed 
in  the  TNF  algorithm,  will  be  easily  obtained  as  functions  of  the  parameters 
of  the  first— stage  linear  filter,  Sq .  Figure  2  shows  a  block  diagram  repre¬ 
sentation  of  the  suggested  scheme. 

An  interesting  broad  class  of  nonlinear  systems  of  the  following  form 
is  considered  here. 


dxt  =  F (x,y , t)dt  +  G(x,y,t)dwt  , 
dyt  =  H (x,y , t)dt  +  a(t)dvt  , 


(15) 


where  w^_,  v  are  mutually  independent  vector  Wiener  processes  of  appropriate 
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Figure  2.  Block  Diagram  Representation  of  the  TNF . 


12 


dimensions,  <j(t)  is  a  matrix  of  compatible  order.  The  functions  F(0,  H(-), 
and  G(*)  can  be  partitioned  as  follows: 

F(x,y,t)  =  fx(t)x  +  f2(x,y,t)  , 

G(x,y,t)  =  gx(t)  +  g2(x,y,t)  ,  (16) 


hx(t)x 

°1  0 

H(x,y,t)  = 

,  o(t)  = 

h2(x,y,t) 

i - 

o 

Q 

K> 

Here  f  h^,  and  g^  are  matrices  of  appropriate  dimensions.  f2>  h2,  and  g2 
are  nonlinear  functions  of  their  arguments,  and  of  compatible  orders. 

The  general  outline  of  the  decomposition  scheme  and  the  various  filtering 
algorithmic  steps  are: 

(1)  The  nonlinear  system  as  in  (15)  can  be  decomposed  into  two  subsys¬ 
tems. 

Subsystem  I 

dxft  =  fl(t)xft  +  S1(t>dwt  » 

(17) 

dyft  =  h1(t)xft  +  a1(t)dvd  , 

where  f^,  h^,  and  g^  are  as  defined  before,  and  w^ ,  are  indepen¬ 
dent  Wiener  processes  of  appropriate  dimensions.  The  subscript  f 
denotes  the  first  subsystem. 
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Subsystem  II 

dxgt  =  f2(x,y,t)dt  +  g2(x,y,t)dwt  , 

2 

dygt  =  h2(x,y,t)dt  +  dvt  , 


(18) 


where  an<^  ^2  are  aS  before  and  wt,  are  again 

independent  Wiener  processes  of  compatible  orders.  The  subscript 
s  denotes  the  second  subsystem. 

(2)  Apply  a  classical  filtering  technique,  i.e.,  the  Kalman-Bucy  algo¬ 
rithm  [6]  to  the  linear  system  in  (17)  as  follows: 


A  A  "1 

dxft  =  f^dt  +  ph*(cr^a*)  dnfc  » 

dpt  =  (flP  +  pfl  +  gl8l  ~  Ph*(a1a*)-1h1p)dt  , 

(19) 

d\  =  dyft  -  Vftdt  , 

xf(0)  =  E(xf (0)) ,  p  (0)  =  cov(xf(0))  , 

where  is  the  estimate,  and  p(t)  is  the  error  covariance  matrix. 
This  will  be  considered  as  the  first  stage  of  the  filtering  algo¬ 
rithm. 


(3)  Find  an  appropriate  bilinear  approximation  model  for  the  system  in 
(18)  of  the  following  form: 
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dxs  ~  f2^x»ut»t^dt  +  g2^ut,t:^dw  ’ 


dys  =  h2^x»ut»t)dt  +  °2dv2  ’ 


(20) 


where 


n 


f2(x,ut,t) 


utixi  +  uc(n+l)  '  A(Vt)x  +  B(uc't)  ' 


n 


h2(x,ut,t)  =  ^  u^x.  +  ut(n+2)  =  C(ut,t)x  +  D(ut,t)  , 


®2(ut,t) 


=  u 


t (n+3 ) 


*  ®o<ut,t^ 


(21) 


Here  the  second  equality  is  used  for  mathematical  convenience.  The 
controls  u  =»  j=l,2, . .  .n+3  are  measurable  with  respect  to  a- 

algebra  {yfg;  s  e  [0,t]},  and  are  chosen  to  minimize  the  following 
global  filtering  criterion 


Q(u)  -  min  E  [  /  (k-k)2  dt]  ,  (22) 

u  0 

where  k  denotes  any  of  the  functions  f2>  or  h2,  while  k  denotes 

the  corresponding  approximation  f2,  g2,  or  h2  in  (21).  Using  the 
property  of  expectation  and  Bayes  formula,  (22)  becomes 


Q (u)  =  min  E(  /  E(k-k)2/y  ;  s  £  t)dt 

u  0  ts 


T 

=  min  E(  f  L(y  ,u)dt)  , 

0  r 


(23) 
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t  ~  2 

where  L(yf,u)  =  E  (k-k)  ,  and  arguments  (x,y,t),  (ut>t)  are  omitted 
for  brevity. 

The  minimization  of  (23)  with  respect  to  u  =  {u^  },  j=l,2,...n+3 

can  be  performed  "locally"  since  L  depends  on  the  cr— algebra  (y^  , 

s  e  [0,t]),  which  is  not  affected  by  the  u,'s. 

3 

Assume  that  the  u^'s  which  minimize  (23)  are  obtained. 

(4)  The  new  equivalent  system  has  the  following  form: 


dxt  =  (A]/Ut,t:)xt  +  +  G1(ut,t)dw  , 

A  A 

dyt  =  +  D1(ut»t))dt  +  CTdv  » 


(24) 


where 


A1(ut,t) 


C1(ut,t) 


G^ (u^ , t) 


[fx(t)  +A(ut,t)],  B1(uj.,t)  =  B(u.,t)  , 

hi(t)  .  0 

,  D  (u  ,t)  =  .  , 

C(ut,t)  1  c  D(ut,t) 


(g1(t)  +  gQ(ut,t)) ,  a  = 


°1  0 


0  O', 


Here  again  the  matrices  are  of  compatible  orders. 

(5)  Again,  with  certain  assumptions  about  (A^,  B^,  C.^  D  ,  G..,  a)  and 

the  distribution  of  the  initial  state  xQ  given  yQ  (see  [2]  and  [3]), 
the  corresponding  conditionally  Gaussian  filter  is  of  the  follow- 
ing  form: 
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dxt  =  ^AlXt  +  Bi^dt  +  S  dv  , 

S  -  (rtC*)(aa*)‘°*5  , 


dv 


(aa*)  (dy ^ 


(ClXt  +  Dl)dt)  ’ 


(25) 


drt  =  (Airt  +  rtAl  +  G1G1  "  ss*)dt  » 

where  A^,  B^,  C^,  D.j^  are  as  in  (24),  and  the  arguments  (U(_,t)  are 
agaxn  omitted  for  brevity.  A  schematic  diagram  representation  of 
the  algorithm  is  given  in  Figure  3.  A  second-order  example  to 
demonstrate  the  algorithm  steps  follows. 


Example  2 

A  second-order,  linear  sonar  target  track  is  considered  here  with  a  noa 
linear  observation.  This  could  represent  active  tracking  of  a  multi-mode 

range  system  or  passive  tracking  with  multi-receiver-transmitter  and  corre- 
lated  time  delay. 

It  is  assumed  that  state  vector  (range  =  x^  range  rate  =  x2)  evolves 
according  to  the  following  stochastic  differential  equations: 


dx  =  Fxdt  +  Gdw1  , 
t  ’ 


dx- 

1 

where  dx  =  j 

1 

1 

i  i 

dx2  I 

-  _J 

,  F  = 

1 

wt  is  a  Wiener  process,  and 


i 

1 

>  G  = 

0 

•a 

a 

- 

j 

is  the  target  maneuver 


1/x. 

time  constant. 


(26) 


Implies  function  or 
fund  tonal  operator 
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Figure  3.  Schematic  Representation  for  the  TNF. 
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The  measurement  equations  are  nonlinear  due  to  target  motion  during 
transmission.  Neglecting  the  nonlinearities  in  the  "velocity  observation,"  it 


is  assumed  that 


dy  =  H (x, t)dt  +  Rdw^ 


(27) 


2  . 

w  is  a 


L  dy2  J  x2  0  ®2_  ' 

Wiener  vector  of  measurement  noises,  m  =  a/c,  a  is  a  constant,  and  c  is  the 


average  speed  of  sound  in  water.  Applying  the  above  algorithm,  the  two 
subsystems  are  given  next. 

Subsystem  I 

Here  only  the  linear,  y ^ 9  observation  is  used.  The  system  and 
observation  equations  are 

4 

dxl  =  x2  dt  » 


(28) 


Then  using  the  Kalman-Bucy  filter  for  (28)  yields 


(29) 
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where  x^  -  E(x^/y2),  i  =  1,2,  the  conditional  expectation,  dv  is  the 
innovation  process,  and  p^,  p^  are  the  standard  time-variant  solutions 
to  components  of  the  Riccati  equation. 

Subsystem  II 

The  best  bilinear  observation  approximation  is  obtained  such  that 

u.,  j  =  1,2,3  are  measurable  with  respect  to  the  a-algebra  (y?  , 

3  zs 

s  e  [ 0 , t ] ) ,  and  minimize  the  following  global  criterion: 

T 

2 

J(u)  -  min  E (  /  (x^  -  (u^  +  u2x2  +  u3)  )dt  ,  (30) 

u.  0 

J 

which  can  be  written  equivalently  using  the  properties  of  expectation  as 

T 

J(u)  =  min  E (  /  L(y,u)dt)  ,  (31) 

u  0 

where  L(y,u)  =  E^x^x2)  -  E^ (R^)  +  Et(N2)  ,  N  -  u^1  +  u^  +  u^, 

R1  "  ^X1X2N‘  ^  re^ers  t0  conditional  expectation  with  respect  to 

the  observation  Y^.  Thus,  performing  the  minimization  with  respect 
to  u  ,  j  -  1,2,3,  the  following  are  obtained  (see  [5]): 

A  A 

U1  =  X2f  ’ 

A  A 

U2  =  Xlf  ’ 

A  A  A 

U3  ~  P  3  Xlf X2f  * 


(32) 
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Then,  the  new  equivalent  system  is 


dx  =  Fx  dt  +  G  dw^  , 


dy  =  C^x  dt  +  D^(x)dt  +  a  dwt  , 


where 


m(l+u1)  u2 


Cl  = 


,  D1(x)  = 


mu. 


0  '! 


_0  a2- 


,  F,G  are  the  same  as  in  (26) 


(33) 


Now,  assume  the  following: 

(i)  if  f  denotes  any  of  the  functions  A,  G,  C^,  D^,  a,  then 
T 

P(  /  |f|  dt  <  »)  -  1  ; 

0 

(ii)  Xq  given  is  conditionally  Gaussian,  Then,  from  (6)  the  corre¬ 
sponding  conditionally  Gaussian  filter  is 


dx.^  =  x2dt  +  (1^(1  +  mu^)  +  mu2r3)dv1  +  — —  dv2  , 

°1  °2 

*  A  1  ~  ^2 

dx2  =  -ax2dt  +  [—  (r3(l  +  mu^)  +  mu2r2)]dv-  +  —  dv2  , 


(34) 


where, 

dv„ 


dy^  -  (1  +  u^mx^  +  mu2X2  +  mu3^t  > 


dv2  =  dy2  "  X2dt  ’ 
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and  m,a  are  defined  as  before. 

The  covariance  equations  are 

p2 

dri  =  {2r3  "  [i2  (ri(1  +  mu2)  +  r3m^2^ 2  +  ~f]>  dt  , 

C1  °2 

2  ,  r2 

dr2  =  {a  -  2a  r2  -  [—  (r3(l  +  mUl)  4-  r2mu2)2  +  -jj}  dt  , 

al  *2 

dr3  =  {r2  "  ar3  ~  (ri(1  +  mui>  +  r3mu2)(r3(1  + 

+  r2mu2>]  -  (r2r3)}dt  ,  (35) 

°2 

A  a  A 

where  u.^,  u2>  u3  are  as  in  (32). 

Notice  that  in  this  case  (34)  and  (35)  are  the  same  as  the  filter 
equations  of  the  modif ied-second-order  truncated  filter  defined 
by  Jazwinski  [6]  because  the  nonlinearity  is  of  second  order,  al¬ 
though  the  approach  is  quite  different. 


Extended  Kalman  Filter 

The  filter  equations  are  [6] 


dXl  X2dt  +  12  +  “*2^1  +  P3mx1Jdv1  +  dv  , 

al  a2 

dx2  =  -axtdt  1(1  +  mx2)p3  +  p2mx1)]dv1  +  dv2 


(36) 


where  m, a  are  as  defined  before,  and 
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dv1  =  (dy1  -  (1  +  mx2)x1dt)  , 


dv2  =  (dy2  -  x2dt)  . 


The  covariance  equations  are 


dpx  =  (2P3  “  (d  +  x2m')pi  +  x1mP3)2  "  “§}  » 

°1  °2 

p2 

dp 2  =  (a2  -  2ap2  -  [(1  +  mx2)p3  +  mx1p2]2  -  -|}  dt 

ai  °2 


dP 3  =  (P2  ~  aP3  “  l~ J  ((1  +  +  x1mP3)((1  +  mx2^ 

al 


p3  +  mp2x  )  -  —  p2p3>  dt  . 
°2 


(37) 
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IV.  SIMULATION  RESULTS 


The  TNE  algorithm  and  the  EKF  algorithm  for  the  previous  example  were 
simulated  by  a  digital  computer.  In  the  simulation,  a  fourth-order,  Runge- 
Kutta  integration  algorithm  was  used  for  all  trajectory  filters  and  differ¬ 
ential  equations  of  both  the  original  system  and  the  error-covariance 
matrices.  Throughout  all  the  simulation  cases,  the  Wiener  processes  w^, 
which  describe  the  excitation  noises,  were  generated  from  pseudo-random 
Gaussian  variables,  v^,  N(0,1).  The  latter  was  generated  by  standard 
(IMSL)  library  subroutine,  and  increments  of  w^  were  approximated  by 
Aw  -  v^At  v^,  where  At  is  the  integration  step  size. 

The  performance  of  the  two  filters  are  compared  on  the  basis  of: 

(1)  The  "mean-square  error"  (m.s.e.)  of  the  filter  output  to  x^ 
such  that 


T 

JTNF  =  E(0i'  "  XTNF(t))  dt)  * 

T 

JEKF  =  ^  ^X^  "  » 


(38) 


and  JJ  gives  the  relative  (percentage)  difference  between  J  , 

I  Nr 

JEKF  SUCh  that 


JJ 


■Win?  x  100 

JEKF 


(39) 


(2)  Root  mean  square  error 
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»  (x<«co -;w(,  »2 


Qj(t)  -  [  i  — I - i-jj-J - L.]1/2  ,  s  -  1,: 

J  i=l  w 


N  (xfi)(t  )  -  x(i)(t  ))2 

VQ.(t)  =  [  z  — *1 - — - ^ 

J  i-1  N 


(40) 


]  '  ,  3  =  1,2 


where  ,  j  “  1,2  are  the  RMS  range  (position)  errors.  VQ.,  j  = 

1,2  are  the  RMS  velocity  (range  rate)  errors,  and  (xf^ (t  ),  xf1^ (t  )) 

th  j  *  j  k 

are  the  j  components  of  the  true  state  and  its  corresponding  TNF, 

EKF  estimates  at  time  t^  on  the  i*"*1  simulation  run,  in  a  series  of 
N  runs.  For  completeness,  some  comments  on  the  filter  initializa¬ 
tion  seem  in  order  here.  Under  actual  operating  conditions  it  is 
extremely  difficult,  and  indeed  rare  due  to  one  reason  or  another, 
to  obtain  reliable  initial  estimates  of  the  state  vector  and  its 
associated  covariance  matrix.  Consequently,  the  following  set  of 
initial  conditions  are  reasonably  chosen.  Throughout,  the  initial 
range  value  is  5000  meters,  while  the  initial  range  rate  value  is 
assumed  constant  and  chosen  from  the  following  set  (50  m/sec,  500 
m/sec,  1000  m/sec) .  The  initial  condition  of  the  estimates  are 
calculated  according  to  the  following  equation: 


x±(0)  =  x±(0)  +  /p±(0)  n±,  i  =  1,2  (4i) 

where  n±  is  a  random  variable.  The  initial  covariance  matrix  is 


P1(0) 

u> 

o 

1 

r  _  _6 
i  10 

t 

102  ' 

i - 

"d 

u> 

o 

V 

P2(°) 

102 

104 

where  the  diagonal  elements  of  p(0)  are  chosen  relatively  large 
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so  that  the  filter  will  "forget"  the  initial  values  as  more  data 
arrived,  and  to  ensure  the  randomness  of  the  initial  estimates. 

In  all  cases,  a  system  noise  of  1%  variance  of  the  initial  state 
values  is  used,  and  different  levels  of  measurement  noise  standard 
deviation  (from  2—20%)  of  the  initial  range,  range  rate  respectively 
are  added.  For  convenience,  the  time  interval  T ,  for  each  run  is 
10  seconds,  and  the  number  of  runs,  N,  for  each  simulation  test 
case  is  20.  Thus,  all  results  have  been  ensemble  averaged  over 
N  =  20  runs. 

The  effect  of  increasing  the  nonlinearity,  (i.e. ,  increases 
in  a),  of  the  system  on  the  rms  error  levels  Q^.(t),  VQ^  (t)  ,  j  = 

1,2,  are  demonstrated  in  Figures  4  and  5  as  compared  to  Figures 
6  and  7,  respectively.  Accordingly,  the  TNF  performance  improved 
substantially,  and  the  rms-error  levels  increased  considerably 
as  compared  to  the  rms-error  levels  of  EKF.  These  comparisons  are 
summarized  by  Table  1. 

Comparison  of  Figures  8  and  9  with  Figures  6  and  7,  respective¬ 
ly  ,  and  Table  2  indicates  that  the  EKF  gains  in  accuracy  relative 
to  the  TNF  as  the  observations  become  more  noisy  (i.e.,  increases 
the  range  measurement  noise  standard  deviation,  a to  20%).  This 
is  due  to  the  fact  that  the  nonlinearity,  (here  in  the  range  measure 
ment) ,  is  masked  by  the  large  measurement  noise. 

From  the  tables  and  figures  mentioned  above,  it  is  seen  that 
in  many  cases  the  TNF  shows  significant  improvement  in  filter 
accuracy  as  compared  to  the  EKF.  For  certain  applications,  the 
complexity  of  the  proposed  algorithm  (TNF)  over  the  EKF  would 
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Figure  5*  Velocity  RMS  Errors  by  TNF,  EKF. 
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Figure  9.  Velocity  RMS  Errors  by  TNF,  EKF. 


Table  1. 


Synopsis  of  the  Percentage  Accuracy  of  TNF  Over  EKF. 


a 

JJX% 

JJ2% 

i 

10.91 

35.57 

1.5 

32.37 

57.14 

2 

78.60 

90.03 

3 

91.42 

98.85 

(cr^  -  2%,  -  10%,  x^(0)  *  5xl0^m,  x^CO)  =  10^  m/sec) 


Table  2 •  The  Effect  of  Measurement  Errors  on  the  Percentage 
Accuracy  of  the  TNF  Over  the  EKF. 


a 

3 

3 

3 

ai 

2% 

10% 

20% 

°  2 

10% 

10% 

10% 

JJX% 

91.42 

11.0 

-  7.23 

JJ2% 

98.03 

17.47 

-16.13 

(x^CO)  =  10^  m/sec,  x^(0)  =  5xl0^m) 
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be  justified  by  the  significant  improvement  in  the  filter  accuracy. 
That  the  EKF  performs  slightly  better  in  the  higher  observation 
noise  case  is  due  partly  at  least  to  the  suboptimal  approximation 


used. 
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V.  CONCLUSION 


A  new  global-filtering  approximation  for  a  certain  class  of  nonlinear 
systems  is  presented.  An  important  practical  feature  of  the  proposed  method 
is  the  method's  independence  of  the  model  smoothness  assumption  which  is 
crucial  to  traditional  techniques.  Furthermore,  a  major  and  equally  important 
byproduct  is  the  generation  of  a  "close"  (in  m.s.e.)  bilinear  model  approxi¬ 
mation  of  the  original  nonlinear  system.  The  assumption  that  Xq  given  y^  is 
conditionally  Gaussian  may  be  satisfied  under  somewhat  realistic  operating 
conditions,  and,  of  course,  it  is  more  general  than  the  traditional  Gaussian 
assumption  of  both  x^.  and  y^. 

The  digital  computer  simulation  demonstrates  the  substantial  filtering 
accuracy  improvement  of  the  TNF  over  the  popular  EKF  in  most  cases. 

In  terms  of  computation  time  the  TNF  takes  roughly  10%  longer  than  the 


EKF. 
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